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ABSTRACT 

Context. The Rayleigh-Taylor instabilities that are generated by the deceleration of a supernova remnant during the ejecta-dominated 
phase are known to produce finger-like structures in the matter distribution that modify the geometry of the remnant. The morphology 
of supernova remnants is also expected to be modified when efficient particle acceleration occurs at their shocks. 
Aims. The impact of the Rayleigh-Taylor instabilities from the ejecta-dominated to the Sedov-Taylor phase is investigated over one 
octant of the supernova remnant. We also study the effect of efficient particle acceleration at the forward shock on the growth of the 
Rayleigh-Taylor instabilities. 

Methods. We modified the Adaptive Mesh Refinement code RAMSES to study with hydrodynamic numerical simulations the evolu- 
tion of supernova remnants in the framework of an expanding reference frame. The adiabatic index of a relativistic gas between the 
forward shock and the contact discontinuity mimics the presence of accelerated particles. 

Results. The great advantage of the super-comoving coordinate system adopted here is that it minimizes numerical diffusion at the 
contact discontinuity, since it is stationary with respect to the grid. We propose an accurate expression for the growth of the Rayleigh- 
Taylor structures that smoothly connects the early growth to the asymptotic self-similar behaviour. 

Conclusions. The development of the Rayleigh-Taylor structures is affected, although not drastically, if the blast wave is dominated 
by cosmic rays. The amount of ejecta that reaches the shocked interstellar medium is smaller in this case. If acceleration were to occur 
at both shocks, the extent of the Rayleigh-Taylor structures would be similar but the reverse shock would be strongly perturbed. 

Key words. ISM: supernova remnants - Physical data and processes: Acceleration of particles. Hydrodynamics 



1. Introduction 



In young supernova remnants (hereafter SNR), the dense shell 
of material ejected by the explosion and decelerating in a rare- 
fied interstellar medium (hereafter ISM) is expected to be sub- 



ject to hydrodynamic instabilities (Gull 1973 1975[ Shirkey 



1 1978 )) of Rayleigh-Taylor (hereafter RT) type. These instabili 
ties modify the morphology of the SNR causing a departure of 
the ejecta from spherical symmetry. They manifest themselves 
as finger-like structures of material protruding from the con- 
tact discontinuity between the two media into the ISM heated 
by the forward shock, as shown by numerical simulations (e.g. 



Chevalier, Blondin & Em meringJ9 92; Dwarkadas & Chevalier| 
1998[|Dwarka das 2000; W ang & Chevalier||2001| l.T)iriing this 



process, the shocked ejecta and the shocked ISM remain two 
distinct fluids. The X-ray observations of Tycho's SNR ( |Warren| 
|et al.||200 5^ exhibit structures that are consistent with these ef- 
fects. Deviations from spherical symmetry may also be caused 
by initial asymmetries in the explosion of the progenitor or lo- 
cal inhomogeneities in the circumstellar medium. In Cassiopeia 
A, the spatial inversion observed by the Chandra observatory of 
the iron and silicon layers (Hugh es et al.|2000| provides a strong 
indication of an asymmetric explosion. An even more radical ex- 
ample of deviation from spherical symmetry is SN 1993 J, a stel- 
lar wind case, where the optical and radio observations can be 



reconciled by assuming a strong asphericity for the pre-existing 
progenitor activity ( Bietenholz et al.|2001 



Two-dimensional hydrodynamic simulations of SNRs can 
take into account the RT instability. Three-dimensional simula- 
tions have shown an enhancement of small-scale structures and 
more severe deformation of the reverse shock suiface Pun &| 
Norman 1996); the perturbation seems to grow faster by 30% 
than in the two-dimensional case ( jKane et al.|2000| ). We chose to 
pursue three-dimensional hydrodynamic numerical simulations, 
focusing on the deviations from spherical symmetry of the SNR 
in the ejecta-dominated phase induced by the RT instabilities. 

Several distinct physical processes occur in young SNRs, 
for instance, the dispersion of synthesized materials in the cir- 
cumstellar medium and the propagation of collisionless shocks. 
Galactic SNRs are also considered to be a strong candidate 
source of galactic cosmic rays up to the knee, i.e. E ~ 3x10'^ eV 
(Lagage & Cesarsky 1983 ), since the rate and energy budget can 
account for the galactic energy density of cosmic rays; however, 
a direct identification of SNRs as sources of cosmic -ray nuclei 
is still lacking. 

In the present paper, we aim to investigate with hydrody- 
namic equations the growth of Rayleigh-Taylor instabilities in 
the context of efficient particle acceleration. We adapt to SNR 
physics the hydrodynamic version of the AMR (Adaptive Mesh 



Refinement) code RAMSES (Teyssier 2002), designed origi 



nally to study the large-scale structure formation in the universe 
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with high spatial resolution. The first application, considered 
here, is to study the effect of the growth of RT instabilities on the 
profile of the hydrodynamic variables by considering the evo- 
lution of a full octant of the SNR, i.e. a larger angular region 
than previously considered (Blondin & Ellison 2001| l. We do 
not introduce any seed perturbation into the 3D radial velocity 
field, so that any departure from spherical symmetry is naturally 
produced by the numerical fluctuations; this is supported by the 
non-linear phase of the instability being insensitive to initial per- 
turbations dChevalier, Blon din & Emmering|1992| l. 

fBlondin & Ellison (2001| described the eflicient acceleration 
of particles by changing the effective adiabatic index through- 
out. A higher compression ratio at the shocks was found to only 
slightly affect the growth of RT instabilities. However, the sug- 
gestion that the reverse shock can efficiently accelerate particles 
continues to be debated ( [Ellison, Decour chelle & Ballet 2005). 
In this paper, we study the impact of cosmic-ray acceleration on 
the instabilities using a simplified description that assumes that 
the adiabatic index is 4/3 in the shocked ISM region, but remains 
5/3 inside the contact discontinuity (surface separating ejecta 
and ISM). This simulates the case of a cosmic ray-dominated 
blast wave and gas-dominated reverse shock. 

To follow the expansion of SNRs, we use synergically two 
numerical approaches: AMR and the Moving Computational 
Grid. As a result, the large-scale turbulence in SNR can be simu- 
lated, allowing for an accurate description of the instabilities. 

The hydrodynamic equations can be formulated with re- 
spect to two distinct classes of coordinate systems: Eulerian, 
i.e., fixed space-coordinates in time, whose major concern is the 
production of numerical diffusion due to the advective terms; 
and Lagrangian, i.e., coordinates comoving with the bulk ffuid, 
which is free in principle from numerical diffusion, but possi- 
ble grid distortions require a rezoning which produces new nu- 
merical diffusion. Therefore, Eulerian coordinate systems are 
generally selected for multidimensional flow simulations. The 
numerical approach of Moving Computational Grid consists of 
using a computational grid comoving, or quasi-comoving, with 
the hydrodynamic flow to minimize the local fluid velocity. The 
idea of a computational grid adapted to follow the global motion 
of the fluid was first proposed in cosmological numerical simu- 
lations by Gnedin ( 19951. However, further analysis ( [Gnedin & 



|Bertschinger|1996 



highlighted problems in the coupling of the 
hydrodynamics so ver with the gravity solver Strong mesh de- 
formations were also found in the gravitational clustering. In the 
approach of Moving Computational Grid, the mesh moves con- 
tinuously and in addition to a full transformation of coordinates 
(position and time), a transformation of hydrodynamic variables 
(density, velocity and pressure) is performed. In contrast, in a 
purely Lagrangian approach the space-coordinates are comov- 
ing with the bulk flow, leaving the hydrodynamic quantities un- 
changed. 

The main disadvantage of simulating a SNR in a Eulerian 
fixed computational grid is that in the bulk flow of the SNR the 
total energy is roughly equal to the kinetic energy. Therefore, the 
thermal energy, computed as the difference between total and ki- 
netic energies, can be very small leading possibly to local nega- 
tive thermodynamic pressure. An algorithm is required to control 
the sign of pressure. 

We apply a combination of the AMR with the Moving 
Computational Grid from the young phase, starting soon after 
the self-similar profile of |Chevalier| (| 1 982 1 983) has been estab- 
lished. The modification introduced enabled us to follow the evo- 
lution of one eighth of the volume of a young SNR, a large volu- 
me with respect to previous 3D simulations (e.g. |Jun & Norman] 



|1996[ ), and for a long interval of time, namely until the transition 
to the Sedov-Taylor phase. We do not provide all the details of 
the RAMSES code, which can be found in Teyssier |(2002) , but 
a description of the modifications introduced here is outlined. 

The plan of the paper is the following. In Sect. |2] we dis- 
cuss the initial conditions adopted for a young SNR, namely the 
mapping of the uni-dimensional solution of the hydrodynamic 
equation over a cartesian grid. In Sect. [3] we present the hydro- 
dynamic equations for the SNR flow, written in both the labora- 
tory frame and the reference frame which is comoving with the 
contact discontinuity. In Sect.|4] we sketch the main steps of the 
implementation of the Godunov method in RAMSES, we dis- 
cuss the modifications and the new variables introduced: the ac- 
tive variable a to change locally the equation of state; the passive 
scalar /, which traces the surface of the contact discontinuity; 
and the ionization age r. In Sect.|5] we present the results of SNR 
simulations: the growth of the RT structures and the elongation 
in time; the influence of a cosmic-ray-dominated blast wave on 
the RT instabilities. In Sect.|6] we discuss the main findings and 
present additional numerical tests. In Sect. [7] we conclude and 
present forthcoming applications. 



2. Initial conditions for self-similar expansion of 
young SNRs 

Supemovae observations have confirmed (for the case of SN 
1987A, 



see 



Arnett 19881 that at early times, typically a few 



hours after the supernova explosion, the matter density of the 
outer ejecta reaches a power-law profile as a function of radius 
given by p = po('"/'"o) "■ A kinematical study of the pre-Sedov- 
Taylor SNR of Kepler (Vink 2008a| l confirmed the power-law 
density profile. The difficulties in attaining a measurement of n 
makes the comparison of theory with observations difficult. The 
value of n depends on the properties of the progenitor star and 
the acceleration of the shock during the explosion. The value 
n - 1, or alternatively an exponential density profile for ejecta 
( jDwarkadas & Chevalier|[T998| l, is commonly used to describe 
Type la supernovae, which are believed to result from the explo- 



sion of a C/O white dwarf accreted in a binary system ( Hoyle & 



Fowler 1960 1. The value n = 9 is believed to describe Type II su- 
pemovae, which result from the core-collapse of a massive star. 
If the density of the ambient medium has a power-law behaviour, 
namely p - pismif / ^0) ^ self-similar solutions can be found for 
the density, velocity, and pressure of the interaction region in 



young SNRs ( [Chevalier 1982 1; from dimensional analysis, the 
time evolution of the r adius of the contact discontinuity has the 
form ( |Chevalier|l982| r(t) oc f-*, where A ^ (n - 3)/(n - s). 

During the expansion, a colder and denser fluid, the shocked 
ejecta, slows down in a hotter and less dense fluid, the shocked 
ISM. Seen in the reference frame of the contact discontinu- 
ity, this deceleration is equivalent to a gravity field pushing the 
shocked ejecta toward the shocked ISM. This instability is at 
the basis of the familiar RT structures (Chevalier, Blondin &| 
Emmering 1992[ l. The structure of the interaction region be- 
tween the ballistically expanding ejecta and a uniform ISM in 
young SNR was first explored numerically by Gull ( 1973 1 975[ l, 
who described how the formation of optical filaments is brought 
back to the RT instability and provided a qualitative explana- 
tion of the amplification of the magnetic field in the interac- 
tion region. Given spherical initial conditions, we investigate the 
growth of RT instabilities in different situations. 

In this paper, a steep power-law density ejecta is assumed to 
expand into a uniform ISM (s = 0). In the innermost part, the 
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common cutoff in radius is adopted to introduce the density 
plateau for the ejecta core 



„ (t A-i Pc(0('-/'-c)"" if r > re 



(1) 



where Pc(0 = g"t"^ r^" oc r , and g" is the normalization con- 
stant in Chevalier (19821. As expected, Pc(0 dilutes from the 
initial high values in the inner core down to a negligible value 
during the Sedov-Taylor phase. 

At a given time to, the physical parameters that fix unequivo- 
cally the profiles of p, u, and P are the total kinetic energy £0 
of the explosion, the total ejecta mass Mgj, the circumstellar 
medium density pisM and the indices of the power-law density 
profile of ejecta and ISM, namely (n, s). The outer radius of the 
inner plateau is given by 



Kr — t , 



(lOn-5 £0 



3 n-3M, 



(2) 



A SNR can be defined to be young when the mass of interstellar 
matter Mism swept up is a small fraction of the mass of ejecta 
Mej given by Mism < A more quantitative criterion is that 
the young SNR phase finishes when the reverse shock radius trs 
has reached the radius of the core of the remnant. 

The self-similar profiles of mass density p, fluid velocity u 
and pressure P for a young SNR dChevaUer 1982 1983[ ) are used 
as initial conditions (see Fig.[T]i. By assuming spherical symme- 
try, the self-similar solution is mapped in a 3D computational 
grid. At the initial time fo = 10 yr, the self-similar profiles are 
assumed to be well-established. It is reasonable to assume that 
such a spherical symmetry has not yet been significantly modi- 
fied by the convective motions grown between the explosion of 
the supernova and the time fo. Given the value of fo, we can fix 
rc(fo) and therefore the density profile is determined down to 
r - Q, where the origin of the explosion is located. We assume 
the reasonable values of SNR parameters £0 = l-6x 10^' erg and 
Mej = 5.0 X Mq, where Mq = 1.989 x lO"'-' g is the solar mass 
and the circumstellar medium density pisM - 0.42 amuxcm"^. 
The only information concerning the progenitor of the SNR is 
therefore contained in the index n and the mass Mgj . For the pre- 
vious set of parameters, we find the initial radii of the contact 
discontinuity, forward shock and reverse shock to be, respec- 
tively, rcD = 0.348 pc, rps = 0.412 pc, and trs = 0.326 pc. 



3. Equation of hydrodynamics in accelerated frame 

3.1. Euler equations in laboratory frame 

In an inertial cartesian (laboratory) reference frame R = {f, r), 
the non-relativistic Euler equations of fluid dynamics, in the ab- 
sence of viscosity and dissipative processes, are written in the 
quasi-conservative form 



d_ 
dt 



( ^ 1 


1 


pUi 

E 




■. a , 


V 



pUj 




( "1 


pUjUj + djjP 







uj{E + P) 














(3) 



where p is the mass density, m,- is the i-th component of fluid ve- 
locity, P is the pressure, E - e + \pu^ is the total energy density, 
and e is the internal energy density. The solution of these equa- 
tions requires the use of an equation of state to find the quantity 
Pip given the quantity e/p ( |Synge||T957| l. We use the common 



short-cut for the adiabatic case, consisting of defining the pa- 
rameter y by means of the equation P - (y - l)e. In the non- 
relativistic case, y can be identified as the ratio of the specific 
heats of the gas (y - 5/3). If the pressure is dominated by the 
relativistic particles, y = 4/3. The active scalar a = (y - 1)"' is 
introduced to treat in a purely hydrodynamic way the eff'ect of ef- 
ficient particle acceleration at the forward shock (see Sect. |5.2| i. 

Equation [5] is written in conservative form, except the so- 
called "quasi-conservative" form of equation for a. If the con- 
servative variable pa is used, higher order schemes of shock 
capturing applied to multicomponent flows can produce oscil- 
lations at the interface of different fluids ( |Shyue|1998 Johnsen 
|& Colonius|2006) . 

For a shock propagating into an ideal gas, the compres- 
sion factor at the s hock can be expressed as p^/p; - {y + 
l)/(y - 1 + 2IM\) ([ZeFd ovich & Raizer|[2002| ), where p2 (pi) 
is the downstream (upstream) density of the shock and Mi = 
[piMj/yiPi]^'^^* is the Mach number in the unshocked medium. 
In the limit of a very strong shock, as in young SNRs, Mi — > 00 
and P2/P1 — > (y + l)/(y - 1), depending only on the thermody- 
namic properties of the medium. 

The passive scalar /, representing the mass fraction of ejecta, 
is initially defined as 



fir) 



1 if r < rcD , 
if r > rcD , 



(4) 



where rcu is the radius of the contact discontinuity. The gradi- 
ent of the function / traces the surface of contact discontinuity, 
namely the finger-like structures arising because of the RT insta- 
bilities, and / is propagated by means of 



dt 



(pf) + ^j{pfuj)^Q. 



(5) 



The X-ray spectra of SNRs usually exhibit significant emis- 
sion lines of heavy elements, such as Fe in Kepler ( jSmith et al. 
1989). Therefore, the computation of the X-ray spectrum of 
SNRs requires an evaluation of the ionization produced by elec- 
trons after the passage of the reverse shock. In this paper, we do 
not compute the ionization state by coupling the hydrodynamics 
equations with the equations for the evolution of ionization and 
recombination (for a ID example of this see Rothenflug et al. 
1994). We follow the ionization age, which allows us to compute 
a posteriori the ionization structure assuming a constant temper- 
ature in the fluid element. 

An additional passive scalar representing the ionization age, 
T, was introduced into the code, defined as 



T(f,r) 



= r ne(f', 



r)dt' , 



(6) 



where fsj, is the time since crossing the forward or reverse shock 
at point r and «e(f, r) is the number density of post-shock elec- 
trons in the heated plasma. As a passive scalar, t does not cause 
any variations in the dynamics of the SNR. The ionization age r 
satisfies the equation 



dt 



(pr) + VjipTUj) = pnMP(r) - BPism) , 



(7) 



where ??(x) is the Heaviside function, f ism is the pressure in 
the unshocked ISM, and we can assume for the pressure P(r) 
in the free expansion region, i.e., r < trs. that P(r) < Pism^ 
therefore Eq. |7] is homogeneous outside the interaction region. 
The constant B is the threshold ratio for the detection algorithm 
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of the interaction region. Since the pressure Pint in the interaction 
region satisfies Pint » f ism, we can choose B = 10^. No failures 
in the detection algorithm have been found for lower values of B 
down to B = 10. 



3.2. Euler equations in comoving reference frame 

Now we consider a non-inertial reference frame R = {f, r). We 
deal with the simplified case of a purely radially contracting ac- 



celerated motion, which was introduced in Martel & Shapiro 
1998| l as supercomoving coordinate system; the more general 



Khokhlov|( 2007|l. The frame R is defined by the transformation 
A: Rh^ R 



r - r/fl , 
dt 

p(r, f) - d^p. 



dt = 



where <p and fi are two scaling parameters. The scale factor 
a{f) must be a non-vanishing twice-differentiable function of 
time only. Once the transformation law of the three fundamental 
quantities, namely length, time and mass, is fixed, the tranfor- 
mation law of all the others can be found. From Eqs.|8]and[9] we 
obtain the transformation law of the velocity vector 



u = «^(u-^r) 



(11) 



where d = da/dt. The velocity u is decomposed into the veloci- 
ty in the inertial frame R and the relative velocity of the bulk 
motion. To establish a computational grid comoving with the 
hydrodynamic flow, we choose the expansion law 



a(t) = aoit/toY 



where we consider the simple case of A being independent 
of time. In the cosmological framework, the transformation 
to supercomoving coordinates maps the uniform and isotropic 
Hubble flow in a matter-dominated universe to a uniform matter 
distribution at rest with constant thermodynamic properties, in 
the adiabatic case, namely y - 5/3. In that case, the aim is to 
focus on deviations from the Hubble flow that drive the structure 
formation in the universe. In the case of self-similar expansion 
of SNR, the ti-ansformation of |MiEteT&~Sfiapfo| ( |Tg98'| l has the 
effect of converting the mesh from Eulerian to quasi-Lagrangian, 
which is exactly comoving with the contact discontinuity while 
the self-similar expansion applies, but does not experience grid 
distortions as in the case of purely Lagrangian mesh. 
From Eqs.[8]and|9] we also have 



P(f,t) 

e(r,f) 



(13) 
(14) 



In 3D, the choice <p - 3, for any value of preserves the in- 
variance of the mass conservation equation. However, as shown 
explicitly in Poludnenko & Khokhlov ( 2007 1, only the choice 



= 3, y8=l 



(15) 



and d - const for a perfect monoatomic gas with y - 5/3 pre- 
serves the invariance of the Euler equation (Eq. |3]l. In this pa- 
per, we report results obtained by fixing (p and /3 in Eq. 15 and 
a{t) = aoit/toY, fixing the value of A according to the physical 
parameters of the specific problem. Hereafter, the choice oq - 1 



is applied without reducing the generality of the treatment. In the 
reference frame R, the Euler equations (Eq. [sjl become 



d_ 
dl 



( ^ 1 


f 


pUi 

E 




. a , 


V 



pUj ^ \ 

pUjUj + 6ijP 
iij{E + P) 





-p{fjUj)g + m{5 - 3r) 



where E - e + \pu^ , H 



case in the pres ence of rotation is treated in Poludn enko &| @ - 



1 d^a 
a df 



a 



2 I da 
~^\df 



i^l and 



, da 



dt^ " 



Here, for a generic function g(t, r) we have 



(8) 








(9) 




(10) 


(Vg), 



1 



-(Vg)r. 
a 



(16) 
(17) 

(18) 
(19) 



In most physical cases, the bulk motion is accelerated; therefore 
in the momentum equation written in the non-inertial reference 
frame a non-inertial force term will appear on the right-hand 



side. In the right-hand side of Eq. 16 both the second line and 
the first term of the third line represent the non-inertial force; the 
second term of the third line, treated as a source term since it acts 
as the gravitational potential term in the cosmological frame- 
work, may create an energy sink, losing the numerical advantage 
of the conservative form (see Sect.[4]for details). 

Since it is an advection equation, Eq.[5]is unchanged by the 
choice of Eq. 15 

d 



j-(pf) + Vj(pfuj) = 0. 



J- J 2) Equation |7] changes as follows 



d 



-Spr) + ^jipTUj) = ph,a-'§(P(f) - BPism). 



(20) 



(21) 



dt 

In this paper, we consider the early phase of SNR expansion and 
the transition to the Sedov-Taylor phase, thus the radiative cool- 
ing is negligible with respect to the total energy content and the 
approximation of adiabatic expansion of a perfect fluid is justi- 
fied. We assume that the fluid mainly consists of fully ionized 
hydrogen, therefore y - 5/3. We notice that for y = 5/3 the last 
term in the third relation of Eq. 16 disappears. If y 5/3, as in 
Sect. |5.2| this term can be regarded as a cooling or heating term, 
representing a microscopic form in which the internal energy of 
a non-monoatomic perfect gas can be deposited in the fluid as 
vibro-rotational degrees of freedom. 

The physical time f is related to the time t for <p - 3 and 
Iby 



to 



1 



al2A- 



1 



1 



(22) 



In the cases considered, A > 4/7 (n > 7 and s > 0), thus the 
factor {2A - 1) ' > for every A. Therefore, the evolution in i 
ranges from f = 0, for t = fo, to f — » {to/a^){2A - 1)"' for f — > oo. 
This squeezing of the total time interval depends on the parame- 
ters of the progenitor of the SNR. We also note that, because of 
Eq. [9] since > 0, the time in R is compressed relative to the 
physical time in R by a factor a^^', reducing the time step of the 
numerical integration. 
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Fig. 1. Angle-averaged radial density profile at f = 50 yr for 
A = 0.57, or («, s) = (7,0), and y = 5/3 with 1, 2, and 3 levels of 
refinement (128^ in yellow, (128-256)^ in green, (128-512)^ in 
black). For comparison, the self-similar profile is shown in red. 
The density profiles are normalized to their values at the forward 
shock. From this diagram, it can be seen that the localization 
of the forward shock and its compression ratio are not strongly 
dependent on the maximal refinement level. 



4. Numerical code 

4.1. Set-up 

The simulations were performed with the hydrodynamics ver- 
sion of the code RAMSES ( Teyssier |2002|. RAMSES was origi- 



nally developed to study, with the AMR technique, structure for- 
mation in the universe at high spatial resolution. The solver is 
based on a second order Godunov scheme. We modified the con- 
stitutive equations of RAMSES as shown in Eqs. 16 and 17 to 



solve the hydrodynamics flow in an accelerated reference frame 
expanding according to Eqs.[8l|9] and 10 The modifications con- 



cern the hydrodynamic solver routines and the external gravity 
routine, because the change of reference frame introduces a non- 
inertial force in the Euler equations, which can be treated as a 
source term in the momentum and energy equation. We also in- 
troduced an additional variable a, treated as a passive scalar in 
Eqs. [3] and 16 to change locally the equation of state and two 



more passive scalars (see Sect. 3.1 for details) 



The simulations were performed over three levels of refine- 
ment, of value 1-1-9, corresponding to a cartesian grid with 
a number of cells 128^', 256-^, and 512^', respectively. The effec- 
tive maximal resolution attained in the interaction region at the 
set-up of the simulation is 0.4/512 pc ~ 8 x 10""^ pc. Our 3D 
numerical simulations were carried out across one octant of a 
sphere. 

The refinement algorithm is applied to the surfaces of the 
contact discontinuity and the two shocks. The surface of the 
contact discontinuity is detected by means of the gradient in the 
tracing function / (see Sect. 3.1 1, rather than the density gradi- 
ent; in this way, the effectiveness of the detection is independent 
of the parameters (n, s). The forward and reverse shocks are de- 
tected by measuring the pressure discontinuity. The advantage 
is that the pressure discontinuity is greater than the discontinu- 
ities in density or radial velocity and independent of the assump- 
tions about the density power-law index in the free expansion 
medium and the equation of state of the gas in the interaction 



region. Figure [T] shows the angle-averaged radial density profile 
at f = 50 yr for A = 0.57, or («, s) = (7,0), and y = 5/3 with 1, 
2, and 3 levels of refinement. The resolution of the RT structures 
is crucially improved by applying the AMR. 

During the expansion of a young SNR, the only dynamically 
interesting region, represented by the interaction region, occu- 
pies a fraction of about 40% (and constant during the self-similar 
regime) of the whole remnant volume. Moreover the Kelvin- 
Helmoltz instabilities triggered by the growth of the RT finger- 
like structures during the "self-similar" phase lead to a strongly 
turbulent phase. Therefore, AMR hydrodynamic simulation that 
attempts to spatially resolve the growth of the RT instabilities 
must refine over an increasingly large relative volume. The con- 
tact discontinuity turns out to be the most computationally ex- 
pensive region of the flow. In contrast, the core region within 
the power-law ejecta, where the value of the uniform density is 
fixed by mass conservation, can be easily handled at the lowest 
possible refinement level. In this context, the AMR technique, 
combined with the approach of Moving Computational Grid, is 
found to be suitable. 

We define the boundary conditions of the conservative vari- 
ables in the corresponding six ghost zones. Since we treat 1/8 of 
the volume of the remnant sphere, in some ghost zones (three) 
the boundary condition is inflow, produced by the advection of 
the interstellar material, and in some other ghost zones (three) 
the boundary condition is reflexive, according to the location of 
the interface of the ghost cell. Changes in the order of the assign- 
ment of the ghost cells significantly change neither the resulting 
intercell flow at the shock nor the physical content of the re- 
sult. In the laboratory rest frame, the ISM is cold and stationary: 
UisM = 0. In the shock frame, this implies that the ISM is ad- 
vected onto the SNR shock surface with a velocity field given 
by UisM - -cfitif (see Eq 
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Therefore, a linear interpolation 
(second order) for the reconstruction of the velocity field at the 
boundaries more clearly describes the inflowing material. 

The time step is controlled using a standard Courant factor 
stability constraint. Based on our expanding coordinate system, 
we also included a time step limitation, where a(t) is not allowed 
to vary by more than 10% over the time step. 

The cartesian coordinate system that we use here is not 
well-adapted to spherically symmetric objects such as SNRs 
and spurious geometrical effects are expected. It has long been 
known that an Eulerian code should take care of its violation of 
Galileian invariance: a turbulent phenomenon observed in a fixed 
computational grid can disappear if it is observed from a compu- 
tational grid moving at constant velocity with respect to the fixed 
grid. The goal of our comoving coordinate system is precisely 
to eliminate this effect. This improvement does, however, have 
the problem, which is common to all shock capturing schemes 
dealing with stationary shocks, of the so-called "odd-even de- 



coupling" and the associated carbuncle phenomenon (Peery & 
Imlay 1988 ; Quirk 1994; Pan doTfi & D'Ambrosio | 20 01|l, namely 
a disturbance in the computation of the flow emerging when the 
shock direction is aligned with one of the grid coordinate direc- 
tions. This effect may be amplified when a flow converges in 
a cell from different directions, as in the case of spherical flow 
on a cartesian mesh. The result can strongly deform the shock 
surface and form ripples or bubble-like regions of local sub- 
density. Local changes to the solver were proposed (Kifonidis 
let al. 2003[ l, based on an algorithm for the detection of the shock. 
Reconstruction of the interface values of the state variables is 
performed, but in the cells detected by that algorithm a more 
diffusive solver is used, such as HLL, HLLC or Lax-Friedrichs 
( |Toro|[T999] ), while in all the remaining regions a conventional 
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Riemann exact solver or an approximate Roe solver can be used. 
Solving this long-standing issue is beyond the scope of this pa- 
per, so we did not implement these possible fixes here. 

4.2. Integration method 

We briefly recall the main steps of the integration method of 
RAMSES ( Toro|| 1999) 1. In the absence of a gravitational poten- 
tial, the three-dimensional hydrodynamic Euler equations can be 
written in the conservative form as (see Eq. |3]l 



5U 

~dt 



+ V-F(U) = 0, 



(23) 



where U has components f/, - {p,pUi,E), F has compo- 
nents Fij(U) = (puj,puiUj + dijP,Uj(E + P)), and both / 
and j have values from 1 to the space dimension c/ = 3. 
The space-time integration is performed by defining a cell 
centered on (xi,yj,Zk) of volume defined by the coordinates 
Vij,k = [(Jc,_i,.^,+ i),(}';_i,}';+i),(Zi_i,Zi+i)] and a time interval 
by At = f"^' - f", which is not constant through the simulation. 
In the present section, the quantity n indicates the n-th time step 
of the integration algorithm. The averaged, cell-centered state is 
defined by 



1 



AxAyAz 



U(x,y,z,f)dxdydz.(24) 



Apart from a smooth flow of the conservative variables (t/)"Jj., 
the discretization of space imposes the solution to a unidimen- 
sional Riemann problem along all the space directions. Inside the 
cell Vij^k, the solution {U)"j^. can be reconstructed with distinct 
algorithms, for instance constant (first order), linear (second or- 
der), to obtain the values at the interfaces {U)"._i {U)"^i 

<%-i,r <^>"#4' ^"'l <^>Ur The inten^en flux 

is then computed by using the values at the interface. The time- 
averaged intercell flux is defined by integrating over the planes 
separating neighboring cells 



AtAxAz 



■ I 2 2 

\.xAz 



i,y,z,f)dtdyAzA25) 



yy+i,z,f)dfdxdz,(26) 



H" 



1 r p4 C'*'' H{x,y,z,^i,t)dtdxAy.{21) 
^xAy j,„ J, , J, , 



UM2 At Ax Ay 



The conservative system can then be written 



(28) 



This integral form remains exact for the corresponding Euler 
system since no approximation has been made. Since the inter- 
cell flux is computed at time t"^^-, the method is second order in 
time. The main changes to RAMSES reported here concern the 
definition of the inflow of interstellar matter into the simulation 
box, as explained above, and the introduction of the source term 



S. RAMSES allows us to define an analytical acceleration vec- 
tor, which is treated as an effective gravity vector We have used 
this option, the grid motion being defined analytically by Eq. 12 
The three independent physical quantities, i.e., length, time, anc 
mass, are rescaled according to Eqs.[8l|9) and 10 All the quanti- 



ties adopted until the end of the present section are computed in 
the R frame. 

In presence of a source term, Eq.[23]becomes 



f .V.F(U).S, 



(29) 



where Cf"jj. denotes a numerical approximation to the cell- 
averaged value of (p,pUi,E) at time t" for the cell {i,j,k), and 
the space indices (i,j,k) and the time index (n) are kept for sim- 
plicity the same in the R frame. The numerical discretization of 
the Euler equations with source terms 5"^^^ is given by: 



UJ,k i.J.k 

At 



pn+l/2 _ p«+l/2 
^M/2,j,k ^i-\l2.j,k ^ 

Ax 

^ij+l/2,k ^i,j-l/2,k 

r}«+l/2 _ rJ«+l/2 

"i,j,k+U2 "i,j,k-l/2 



Az 



_ r.n+1/2 
- ^ ij,k 



(30) 



across cell 

j+l/2,k' 1/2 



The time centered fluxes F"^}j^.,, G"^'{^, 
interfaces are computed using a second-order Godunov method 
According to Eq. [16] the gravitational source terms are given by 



,n+l/2 
'i,j,k 





Pa.(f^)a.+p-i(r^)-i 






^i,j,k 



{2 -31a)} 



(31) 



where Q is defined in Eq. 17 and e is the internal energy density. 
The second term on the right-hand side of Eq. 3 1 is equal to zero 
if y = 5/3 or, equivalently, a - 3/2. 



Although the integration of the source term in Eq, 

is of second order in time, higher order methods, such "as 
Rosenbrock ( Kaps & Rentrop|1979| l, can be chosen. 
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5. Application to SNR 

5.1. Growth of Rayleigh-Taylor structures 

Many developed codes have satisfactorily tackled the growth of 
RT instabilities as a mandatory validity test. In contrast to pre- 
vious treatments, we consider the evolution of a full octant of 
the SNR. Thus, within a larger SNR volume, a large number of 
RT-instabilitiy lengthscales, spread over a wider range, sum to- 
gether. In other words, the growth of RT instabilities is more ac- 
curately described not only with high spatial resolution, but also 
with volumes comparable to the global SNR volume. In particu- 
lar, since the initial departures from the spherical symmetry are 
sampled over a larger surface, the initial phase of the RT growth 
is statistically more accurately described (see Sect. 6.3 1. 



Even in the case of a stationary advective shock, a small vol- 
ume has been dealt with, in the SNR literature, in comparison 
with the volume of the interaction region. The description here 



Fraschetti R: Growth of the 3D Rayleigh-Taylor instabiUty in Supernova Remnants using expanding reference frames 



7 



1.6 
1.4 
1.2 
1 

^0.8 



numerical 
TM 

Ejecta-Dominated phase 




Fig. 3. Temporal evolution of the forward shock and reverse 
shock radii for A = 0.57, or (n, s) = (7, 0), and r = 5/3. Our 3D 
numerical solution (solid) matches the ID semi-analytical solu- 
tion ( Truelove & McKee,1999| ), dashed, and, in the asymptotic 
regime for f ^ 0, matches the self-similar solution in the ejecta- 
dominated phase ( ChevaLier|1982| l, dotted. See Sect. 5.1 for the 
units rch and tch {rch ~ 7.8 pc and tch ~ 1950 yr for the cho- 
sen parameters). The spherical symmetry of the reverse shock is 
strongly modified by the RT instability; here the innermost value 
of the radius of the reverse shock has been chosen (cf. Sect. |5.1| l. 
The ejecta-dominated solution (dotted) is prolonged indefinitely 
to emphasize the inward motion of the reverse shock. 



provides a more comprehensive overview of the dynamics of RT 
instabilities in the SNR. 

The simulations presented here cover three phases of the adi- 
abatic SNR expansion (see Fig. |2|: the self-similar expansion 
( Chevalier||1982 1 with an exponent A depending on the proper- 
ties of the ejecta and the ISM, i.e. on («, s); the second phase, 
not self-similar, representing the transition to the Sedov-Taylor 
regime (Mism ~ M^j); and the third phase (Sedov-Taylor), with 
a blast wave radius Rj, that expands as Ri, oc f^/^, for s = 0, 
and Mism M^y The transition between the two self-similar 



phases is shown by means of various snapshots of the angle- 
averaged density in Fig. [2] In this case, we assume that A - 0.57, 
or (n, s) = (7, 0), and y = 5/3 in the whole simulation box. The 
inward motion of the reverse shock and the progressive estab- 
lishment of the Sedov-Taylor profile is manifest. We extended 
our computation until values of time t when the inward motion 
of trs has reached the inner core, and show that during the tran- 
sition to the Sedov-Taylor phase, the angle-averaged profile is 
not significantly altered by the development of the RT instabil- 
ity. The ID Sedov-Taylor density profile is superimposed in the 
last three panels of Fig. |2j by assuming that the Sedov-Taylor 
phase is fully established at f = 5 x fsT, where fsx - 0.732 x tch 
( [Truelove & McKee|1999| l. Compared with the original version 
of the code in the case of a spherical Sedov blast wave (Teyssier' 
||2002 ]l, we note the improvement in the resolution of the com- 
pression factor; in that previous work relatively small values of 
final output times had been considered, the aim being to test the 
code. 

In Fig. |3] the outermost value of the position of the for- 
ward shock and the innermost value of the position of the re- 
verse shock in 3D are compared with the semi-analytical ID 
solution of |Truelove & McKee| ( |1999| l. The units used here cor- 



respond to re, (pc) = 3.O7(Mej/M0)i/ViSM//iH) '^^ and t^h (yr) 

= 423£5;^2(Mej/Mo)5/6(pISM///H)-'^^ where A/H = 2.34x10-^4^ 
is the mean mass per hydrogen nucleus assuming cosmic abun- 
dances. In particular, the inward motion of the reverse shock is 
compared with the solution of Truelove & McKee ( 1999[ l; in the 
case (n, s) - (7, 0), we find agreement within 7% between our 
3D simulation and the ID solution of Truelove & McKee (1999). 
The radius of the reverse shock is strongly modified by the RT 
instability where the innermost value is chosen. In contrast, the 
radius of the forward shock is unequivocally defined, since it is 
not perturbed by the RT instabilities for the physically relevant 
values of (n, s) (see also Fig.[5]l. 

In Fig. |4] (upper panel), the mass fraction /(r), defined in 
Eq. |4] is shown in a planar slice parallel to the coordinate plane 
XY at time t - 100 yr (left) and t - 500 yr (right); in this 
case, A = 0.57, or (n, s) = (7,0) and y = 5/3. The mixture of 
the two media, shocked ejecta and shocked ISM, as well as the 
deformation of the contact discontinuity surface appear at inter- 
mediate values of /(r). In the middle panel, the matter density is 
shown with the same spatial extension and times as in the upper 
panel. In the lower panel, an estimate of the frequency emis- 
sivity 7b produced by bremsstrahlung in purely hydrogen gas is 
shown, integrated along the line of sight towards the reader. Only 
the interaction region where the temperature is high contributes 
to the emissivity. The emissivity /b is computed based on the 
assumption that the plasma is optically thin (bremsstrahlung). 
Therefore, the emissivity is given by A = 2 x 10^^"^ Z^nJ'TI^^, 
where A^e is the electron number density and Tf. is the elec- 
tron temperature given by the perfect gas equation of state such 
that Te = iJ.muP/{pk), where ju is the mean molecular weight, 
me = 1.67 X lO^^^^g is the proton mass, and k = 1.3807 x 10"'^ 
erg/K is the Boltzmann constant. For simplicity, in this case we 
chose yu - 1/2, corresponding to the case of fully ionized hy- 
drogen. We implicitly assumed that Tg = Ti, where Ti is the ion 
temperature, while at an early stage the SNR is not yet com- 
pletely thermalized. A more detailed study should account for 
inclusion of the evolution of T^jTi, but this indicative diagram 
provides nevertheless a starting point. 

In Fig. |5] a 3D large view image of the density field in the 
full simulation box is shown. This picture shows that while the 
sphericity of the forward shock surface is maintained since the 
RT instabilities do not reach the shock, the reverse shock surface 
is strongly modified. The large relative volume of the refinement 
region can be clearly seen. 

In Fig.|6j the streamlines of the velocity field in the reference 
frame of the contact discontinuity in distinct snapshots show the 
development of the chaotic flow in the interaction region due to 
the convective motions. As the time proceeds, the central den- 
sity decreases and the relative size of the region of turbulence 
increases. 

For a general acceleration field g, the initial growth of small 
amplitude perturbations at the accelerated contact discontinuity 
between two incompressible fluids is exponential with a growth 
rate cr. The quantity cr is related to the wavenumber k by the re- 
lation cr= -iJkgA, where A = (pi-p2)/(pi+P2) andpi,p2 arethe 
densities of the two adjacent media at the contact discontinuity. 
At later times, the growth of the ripples at the contact discon- 
tinuity enters the so-called "self-similar phase"; this phase was 
first quantitatively analyzed by Fermi & von Neumann (1953|. 
The self-similar growth of the RT instable structures is governed 
by the equation 



dH(t) 
dt 



2(,5Agfl^H^I^{t) ., 



(32) 
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Fig. 5. The three-dimensional density in the full simulation box 
of side L - 4.2 pc is shown at f = 500 yr for A - 0.57, or 
(n, s) - (7, 0) and y - 5/3. The black line indicates the height of 
the planar slice at z = 0.1 x L used in Fig.|4j upper and middle 
panels. 



where the constant <5 is a dimensionless growth parameter de- 
pending on the dimensionality and the geometry of the system, 
and g is the acceleration. If ^ = const, the asymptotic quadratic 
law in time is found to be Hi{t) = 6Ag(t - fo)^. A second order 
polynomial fit was found to reproduce quite satisfactorily the 
mixing of ejecta and ISM in a young SNR ( Dwarkadas]|2000[ l. 
However, in view of the self-similar deceleration of the con- 
tact discontinuity of an SNR, it is reasonable to assume a time- 
dependent g. We propose that a more accurate result can be 
found by assuming that the acceleration g i s also self-similar. 



namely g{t) ~ ^. The integration of Eq. I32jgives readily 



Hit) A5A{A - 1) 



'"CD(fo) 



A 



A/2 



(33) 



where rcoito) is the initial radius of the contact discontinuity. 
The solution in Eq. [33] reproduces our simulation (see Fig. |7]i. 
The comparison in Fig. [7] shows a disagreement of a factor of 
more than 2 with the second order polynomial fit at early times, 
because the deceleration of the contact discontinuity has not 
been properly taken into account. 

The limited spatial resolution attainable at the contact dis- 
continuity with a cartesian grid makes it hard to verify early time 
exponential growth ( Fermi & von Neumann|1953 i, prior to the 
self-similar regime. 

5.2. Non-thermal particles in the shocked ISM 

Multiwavelength surveys have inferred that cosmic -ray acceler- 
ation can be efficient in young SNRs (see Vink 2008b and ref- 
erences therein). The idea at the basis of particle acceleration 
at SNR shocks is the Fermi first-order mechanism of accelera- 
tion, namely that particles accelerate by repeatedly crossing the 
shock and colliding at every passage with magnetic turbulences. 
Therefore, beside other parameters, the efficiency of acceleration 
depends on the intensity and the spatial configuration of the mag- 



exact solution 

polynomial ■ 



Fig. 7. The protruding RT structures measured in our simulations 
at different times are quite well reproduced by the exact solution 
H(t) with A = 0.57 (see Eq. 33 1. The edge of the RT structures is 
defined by f(r) ~ 0.01. The early-time exponential growth can 
hardly be verified with the available spatial resolution. A second 
order polynomial fit is also shown for comparison. 



netic field. The X-ray detections of narrow synchrotron emitting 
filaments at SNR shocks indicate magnetic field intensities up 
to 1 mGauss (Bam ba et al. 2003 ; Lazendic et al. 2003, \^nk&l 
Laming 2003 ). Magnetic field amplification at shocks has been 
interpreted as the result of turbulence in the postshock material 
(Giacalone & Jokipii 20 07)1 or cosmic-ray induced streaming in- 
stabiHty (Bell 2004; Amato & Blasi 2009). 

The SNR hydrodynamics can be modified by efficient par- 
ticle acceleration at the forward shock ( |Decourchelle, EIlison| 
& Ballet|2000 ). Therefore, a complete numerical description of 
the SNR expansion should take into account the feedback of the 
non-thermal population of particles on the dynamics of the SNR, 
i.e., the system of Euler equations (Eqs.[3]l should be coupled to 
the equation for the cosmic-ray phase-space distribution func- 
tion. 

[Blondin & Ellison) ( [2001) 1 investigated this scenario by 
changing the adiabatic index -y globally. They found that the de- 
velopment of the RT instability was relatively unaffected and the 
main effect was that, because the compression ratio increases 
and the shocked region shrinks when y decreases, the RT fin- 
gers approached the forward shock and could even reach it for 
-y - 1.1. We adopt a similar approach but investigate when par- 
ticle acceleration is efficient at the forward shock only. We call 
this the hybrid case in the remainder of the paper A local equa- 
tion of state of the form P - (yir) - l)e is assumed. The shocked 
ISM is initialized as a relativistic gas (dominated by cosmic rays) 
and the shocked ejecta is initialized as a non-relativistic gas such 
that 



5/3 if r < rcD 
4/3 ifr>rcD, 



(34) 



the corresponding a then being propagated according to Eq. 16 



The initial adiabatic y distribution in Eq.)34]should provide a 
more realistic representation of forward shock accelerating par- 
ticles than a 7 = 4/3 uniform distribution across the whole in- 
teraction region. There has been relatively little observational or 
theoretical evidence that particle acceleration is very efficient at 
the reverse shock ( [Ellison, Decourchelle & Ballet|2005) ). An in- 
teresting result of that simulation is that the density behind the 
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Fig. 8. Angle-averaged radial density profile shown at distinct 
ages for A - 0.57, or (n, s) = (7, 0) in the hybrid case in black, is 
compared with the corresponding initial profile, in red. 



reverse shock is lower than behind the forward shock in the self- 
similar solution, at least for « = 7 (Fig.[8]as opposed to Fig.[T]i. 
This may hinder the development of the RT fingers. 



6. Discussion 

6.1. Intrinsic issues in fhe stationary advective sliock 
problem 

In simulating stationary advective shocks, numerical instabilities 
related to the "carbuncle phenomenon" can be found by using 
the Roe solver. The carbuncle phenomenon was first observed 
by Peery & Imlay ( 1988, ) in blunt body simulations using Roe's 
method and numerically studied in detail by Quirk| ( 1994 1, and 
has remained a crucial numerical issue CPandolfi & D'Ambrosio' 
Xu & Li 200 l| l until recently (Nishika wa & Kitamura 



2001 



2008 L The carbuncle phenomenon is considered to be a spu- 
rious numerical effect affecting hydrodynamic simulations when 
the intercell flow is aligned with one of the coordinate direc- 
tions of the computational grid; it shows up when a supersonic 
shock, or more generally a discontinuity in the flow variables, 
is comoving with the computational grid. The Riemann solver 
produces unphysical ripples of the shock and in the post-shock 
region, which at late times can modify the structure of the shock 
surface itself. The numerical methods that are likely to produce 
carbuncle phenomenon have not yet been unequivocally identi- 



fied, an attempt to classify them being Dumbser, Moschetta & 
|Gressier| ( |2003l ). 

The so-called "H-correction" ( |Stone et al.|2008| ), which con- 
sists of introducing an anisotropic dissipative flux in the direc- 
tion perpendicular to the direction of propagation of the shock, 
cannot be applied here. We did not add a diffusive term in the 
intercell flux, which is usually adopted to avoid the growth of 
numerical instabilities. The addition of artificial viscosity could 
damp physical effects emerging during the development of the 
RT instabilities, which we are interested in studying in detail. 
For the same reason, we did not use more diffusive solvers such 
as HLL or HLLC or Lax-Friedrichs (Toro'1999 ). 

In the case of supernova core-collapse simulations, a possi- 
ble solution has been to compute the intercell flux in the zones of 
strong shock with an Einfeldt solver ( |Quirk|1994| |1997) l, while 



all the other grids in the computational box are treated with a less 
diffusive Riemann solver ( Kifonidis et al.|[2003 i We preferred 
not to introduce a solver-switching algorithm to avoid introduc- 
ing other spurious effects possibly arising from the matching of 
the intercell fluxes. 

To avoid the numerical instability occurring at the stationary 
shock, we used the exact solver The failure of the Roe solver 
manifests itself in local fluctuations of density in regions behind 
the shock and strong deformations of the shock surface, which 
grow as the non-linear regime of RT instabilities is reached; 
consequent unphysical effects due to velocity shear, resembling 
Kelvin-Helmholtz instabilities, develop at the same time. Both 
effects are absent if the shock is not stationary with respect to 
the computational grid and depend on the chosen simulation pa- 
rameters, for instance CFL number or interpolation scheme of 
conservative variables; thus they must be attributed to the choice 
of the solver The choice of the exact solver did not significantly 
increase the computational time. 



6.2. Quasi-comoving grid 

A different solution to the carbuncle problem could involve 
defining a computational grid that drifts with respect to the SNR, 
slowly enough to follow the expansion of the remnant for a long 
time and to avoid boundary reflections possibly occurring with a 
remnant overtaking the grid. However, this is an even less viable 
solution since, as shown in the following, the angle-averaged 
density profiles strongly deviate from the expected compression 
ratio at the two shocks. In contrast, as shown earlier, in the case 
of an exactly comoving computational grid, the post-shock rip- 
ples average to zero. 

Independently of the initial density profile, i.e., of the value 
of A, the self-similar regime slows down towards the Sedov- 
Taylor solution, which has a self-similar exponent equal to 
2/5 < A. Therefore, the contact discontinuity in the computa- 
tional grid will progressively acquire a drift velocity depending 
on the physical parameters of the remnant. To test the numerical 
stability of the shock, we introduced a "quasi-comoving" grid, 
whose small drift with respect to the contact discontinuity ve- 
locity is parametrized by e <k 1: a(t) - a(){t I t^Y^'^ . 

The quasi-comoving grid, even with the Roe solver, smears 
out the ripples produced behind the forward shock surface by the 
stationarity of the computational grid, therefore locally is odd- 
even instability-free (see Fig. |9|l. However, the angle-averaged 
density profile, even at the first time-steps, is affected by an er- 
ror in the compression factor at the two shocks (see Fig. 
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up- 
per panel), the position of the shock being unchanged. The re- 
distribution of the mass through the two shocks in the presence 
of a grid drift (over-compression at reverse shock and under- 
compression at forward shock) is to be attributed to the mapping 
of a spherically symmetric object on a cartesian grid. It man- 
ifests itself during the transition to the Sedov-Taylor phase in 
terms of an error of 5% in the resolution of the compression at 
the shock, and should be considered an intrinsic limit of the ap- 
proach used here. This error does not depend on the number of 
levels of the AMR (see Fig. 10 lower panel). As a consequence, 
even if a larger number of computational cells is considered, the 
angle-averaged value at the shock does not fulfill the value of 
the compression predicted by the Rankine-Hugoniot conditions. 
In contrast, for an exactly comoving grid, even with the Roe 
solver, strong deformations of the shock surface are triggered 
by the numerical instability but the compression at the shocks 
is preserved. The same difference of results between the exactly 
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Fig. 10. Angle-averaged radial density profiles between the re- 
verse shock and the forward shock at time t - 12 yr = 0.0063 
tch, with A = 0.57, or {n,s) = (7,0), and r = 5/3. The AMR 
is performed over three levels of refinement. The radial coordi- 
nate is normalized to the forward shock position corresponding 
to e = 0; the density is normalized to the value at the shock. 
Upper panel. The profile in the case of the exactly comoving grid 
(s = 0) is compared with the profiles having a quasi-comoving 
grid s = (0.0286,0.0536,0.129). The initial profile, at to = 10 
yr, is also depicted for comparison. It is evident that even small 
values of s produce oscillations and an error of the order of 10% 
in the compression factor for e = 0.129. Outside the interaction 
region, the profiles are superimposed. Lower panel. The curves 
without AMR and with AMR are compared with the initial pro- 
file at f() = 10 yr In this case, s - 0.129. The independence of 
the error in the compression factor at the two shocks from the 
space resolution is clear. 



comoving and quasi-comoving has been found with the exact 
Riemann solver, shown in Fig. 10 



As a consequence, during the transition to the Sedov-Taylor 
solution, the generation of a numerical error must be expected 
in the compression at the shock. A possible solution is the intro- 
duction of a grid velocity given by a{t) - a^it 1 
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Fig. 11. Angle-averaged radial profile of the mass fraction of the 
ejecta / shown at distinct ages for A - 0.57, or (n, s) - (7, 0); 
the cases with uniform y - 5 13>, y - A 13 and the hybrid case 
are compared. The radius is in units of the corresponding radius 
of the contact discontinuity. The vertical lines on the right of 
the panels indicate the corresponding positions of the forward 
shock, which at this age still meets the self-similar expansion 
(see Fig.[3]for the uniform case). The vertical line on the left of 
the panels indicates the position of the reverse shock. 



6.3. Cosmic-ray-dominated blast wave 

For i = 0, a convenient case to consider when analyzing the 
growth of RT structures is n = 7, for which the width of the inter- 
action region is large, thus the elongation of the instabilities can 
be studied for a longer time and the forward shock remains unaf- 
fected. The effect of changing the equation of state on the growth 
of the instabilities is represented in Fig. 



1 1 It shows how far the 



RT fingers reach during their growth (100 yr) and when they are 
fully developed (500 yr). At 500 yr, the fraction of ejecta mate- 
rial transported well into the shocked ISM (at r/rco = 1-05, say) 
is lower in the hybrid case. Nevertheless, the outermost ejecta 
reach nearly as far, and very close to the forward shock. 

In the hybrid case, the position of the reverse shock is un- 
affected by a relativistic gas at the forward shock, and the de- 
parture from the self-similar position is less than 2%. If the 
relativistic gas is uniformly distributed across the whole inter- 
action region, the departure of the reverse shock from the self- 
similar position is of the order of 6%. Similar deviations have 
been found by |Blondin & Ellison| ( |200l] ), who considered how- 
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Fig. 12. The protrusion of RT structures measured in our simu- 
lations in units of the corresponding self-similar position of the 
contact discontinuity, whose surface is directly modified by the 
instabilities; the uniform 7 = 5/3 and the hybrid cases are com- 
pared (/(r) ~ 0.1). The cosmic rays affect the elongation of the 
instabilities from the very beginning by a factor of the order of 
20%. 



ever a relatively small angular region in the cases (n, s) - (7, 2) 
or (n,s) = (12,0). The bumps in f(r) found at r/rco ~ 1-05 
originate in the RT structures themselves, independent of the y 
distribution. In the two cases y = 5/3 and y = 4/3, the width 
of the shocked ISM region (rp^/rcD - 1) shrinks from 0.18 to 
0.10. Comparable shrinking has been found in the angle-average 
analysis of X-ray observations of Tycho's SNR (Warren et al. 
|2005] l 



As shown in Fig. 12 the elongation of the RT structures in 
the hybrid case decelerates relative to the non-relativistic case. 
The deceleration is caused by the higher density gradient ahead 
of the contact discontinuity (cf. Fig. [8] and Fig. TJ. This can be 
easily understood because the RT growth rate cr depends through 
A on the average ratio of the two densities in the interaction 
region, i.e., Pej/piSM (cf. Fig.lSll, thus lower density ratios pro- 
duce lower growth rates (see Fig. [T2| . The higher compression 
in the ISM shocked region in the hybrid case produces the slow- 
ing down s hown in Figs. [TT| and \T2\ When y = 4/3 on both 
sides as in Blondin & Ellison ( 2001| , Pej/piSM remains larger 
than 1 and the Rl' hngers reach further than in the hybrid case at 
f(r) < 0.05. 

The region which drives the instability at very early time is 
that part of the density profile where density decreases outwards 
(see Fig.[T]i. This is entirely on the ejecta side and therefore unaf- 



fected by changing y in the shocked ISM (see Fig. 11 1. In other 
words, at a very early stage, when the instabilities grow expo- 
nentially, the process occurs locally in the outer ejecta and the 
cosmic rays cannot affect it. However, particle acceleration will 
be able to affect RT structures as soon as they reach the shocked 
ISM, since the density gradient of the shocked ISM is higher for 
higher compression ratios. For the reasons explained above, a 
phenomenological equation of the type of Eq.[32]can only qual- 
itatively reproduce the global behaviour. We could not find a 
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more accurate description in this case. A comparison of Fig 
with Fig. 7 of |Blondin & Ellison ' (2001) might indicate that the 
large fluctuations they observed in the late-time radial extent of 
RT structures are not an intrinsic property of the instability but a 
result of the too small physical size of the region considered. 



7. Conclusions 

We have adapted the AMR code RAMSES, which is based on a 
second order Godunov method in an expanding frame called "su- 
percomoving coordinates", to follow the evolution of SNRs. In 
this approach, not only the space-time variables have been mod- 
ified but also the hydrodynamics variables. The comoving coor- 
dinate system allows an eighth of the total volume of the SNR 
to be described, which is much larger than previously consid- 
ered. A longer time interval can be investigated, of the order of 
thousands of years, until the transition to the asymptotic Sedov- 
Taylor phase. Such a large volume allows the convective instabil- 
ities to be modeled more accurately since it considers not only 
the shortest wavelengths, which have the greatest growth rate, 
but also the longer wavelengths, which grow more slowly but 
make a significant contribution to the morphology of the SNR 
over time-scales of thousands of years. Our larger spatial sam- 
pling allows a statistically more accurate description of the in- 
stabilities. 

The great advantage of the method adopted here is that it 
minimizes the velocity of the fluid relative to the grid, and the 
numerical diffusion at the contact discontinuity, where the in- 
stability should develop. In the comoving reference frame, the 
contact discontinuity, in the absence of distortions due to con- 
vective instabilities and before the transition to the Sedov-Taylor 
phase, would be exactly stationary in time. The analytical non- 
inertial term resulting from this coordinate transformation is 
strictly equivalent to a gravitational acceleration: the Euler equa- 
tions have therefore been integrated with the corresponding ad- 
ditional source terms. 

The elongation of the Rayleigh-Taylor structures slowly 
reaches the asymptotic behaviour r', by direct solution of the 
self-similar theory according to the assumption of a self-similar 
acceleration instead of a constant acceleration. 

We have presented a simple way to numerically investigate 
the effect of efficient particle acceleration at the forward shock, 
approximating the shocked ISM as a relativistic gas and the 
shocked ejecta as a non-relativistic gas, with a different adiabatic 
exponent in each fluid. The density behind the forward shock 
may be higher than at the reverse shock in that case. A decelera- 
tion in the protruding of RT structures is caused by the higher 
compression of the shocked ISM. The conclusion of Blondin| 
& Ellison (2001 1 that the RT fingers can travel very close to 



the shock in the presence of accelerated particles remains valid. 
The elongation of the instabilities has here been more precisely 
quantified. We propose that this will allow us to understand why 
ejecta are found very close to the blast wave in the remnant of 
SN 1006 ( |Cassam-Chenai et al.|[2008| . The set-up of the code 
presented here will allow future studies of the back-reaction of 
particle acceleration on the SNR evolution. 
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Fig. 2. Angle-averaged radial density profile at various ages for A - 0.57, or («, s) = (7,0), and -y - 5/3 with 1 and 3 levels of 
refinement (128"^ in yellow, 128^ - 512^ in black). For comparison, the red line shows the self-similar profiles normalized to their 
forward shock values. Chevalier profiles for the first five panels and Sedov for the last four panels. The departure from the ejecta- 
dominated phase, occurring when the reverse shock encounters the plateau outer radius r^, and the transition to Sedov-Taylor phase 
are shown. The inward motion of the reverse shock, initially only in the comoving frame of the contact discontinuity and later on in 
the laboratory frame as well, is also depicted. 
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Fig. 4. Snapshots at two distinct ages (f = 100 yr at left; f = 500 yr at right) for A - 0.57, or (n, s) - (7,0) and y = 5/3. In the upper 
panel, the mass fraction fir), defined in Eq. |4] is shown in a comoving slice plane parallel to the coordinate plane XY at height 
z = 0.1 X a(t). The mixture of the two media, namely ejecta shocked and ISM shocked, as well as the deformation of the contact 
discontinuity surface is shown. In the middle panel, matter density is shown in the same plane as in the upper panel. Density values 
are coded according to the color bars given for each frame. Note the change in the radial scale. In the bottom panel, frequency 
emissivity Ji, due to bremsstrahlung integrated along the line of sight towards the reader is shown, assuming a mean molecular 
weight = 1/2, corresponding to the case of fuUy ionized hydrogen. 
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Fig. 6. Snapshots showing streamlines of velocity field superposed to matter density for A - 0.57, or (n, s) - (7,0) and y = 5/3 in 
the frame of contact discontinuity. The planar slice is in all panels z-Q. From left to right, in above panel {t - 100 yr) the inflow 
motion of the ISM and the outflow motion of the freely-expanding ejecta are clearly evident. In the other panels (f = 500 yr, t - 1000 
yr, t - 2000 yr), the turbulent flow motion due to RT instabilities is shown. Note that the scale in color bar is log in amuxcm"^ in 
the first panel and linear in all the other panels. 
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Fig. 9. Matter density in a comoving planar slice parallel to the coordinate plane XY at height z = 0.1 x a(t) (t - 100 yr at left; 
t = 500 yr at right) for A - 0.57, or « = 7, i = 0, and y = 5/3. Density values are coded according to the color bars given for 
each frame. Note the change in the radial scale. In this case, e = 0.129. Since A' > A, the computational grid expands faster than 
the SNR and the volume fraction occupied in the simulation box by the SNR decreases. It is evident that even a small drift of the 
computational grid from the shock velocity erases the effects of the numerical instability possibly emerging behind the shocks. 



